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Introduction 


Breast  cancer  is  the  second  most  common  cancer  type  affecting  American  women.  The 
disease  is  also  the  second  leading  cause  of  cancer-related  death  for  American  women, 
which  was  predicted  to  kill  40,460  women  in  2007.  1  Presently  there  is  no  effective  way 
of  preventing  the  disease.  However,  the  detection  of  the  cancer  at  its  early  stage  has  been 
found  to  significantly  improve  the  survival  rates  2  5.  Moreover,  there  are  more  treatment 
options  for  a  patient  with  earlier  stage  cancer  than  whose  with  a  late  stage  of  cancer  6"8. 

Film- screen  X-ray  mammography  is  still  the  only  FDA  approved  screening  tool  aiming  at 
early  detection  of  the  breast  cancer.  While  it  has  been  proven  to  be  effective,  it  is  not 
omnipotent  in  its  detection  sensitivity  of  breast  lesions  due  to  several  limitations  such  as 
two-dimensional  (2D)  projection  data  acquisition  and  restricted  range  of  linear  optical 
response  of  the  detector.  For  women  with  dense  breasts,  the  sensitivity  is  lower  since  in 
their  mammograms  the  dense  appearance  of  the  breast  tissue  is  more  likely  to  obscure 
any  abnormalities  and  makes  the  detection  of  breast  cancer  even  more  challenging  9. 

With  the  development  of  flat-panel  detectors  in  recent  years,  the  dedicated  cone-beam  CT 
technology  became  feasible.  By  tomographic  reconstruction,  the  issue  of  tissue 
overlapping  is  solved.  The  technology  is  potentially  advantageous  for  women  with  dense 
breasts.  A  research  group  in  University  of  California  Davis  is  currently  conducting 
clinical  trials.10  The  reconstructed  volumes  of  human  subjects  provide  exciting  new  3D 
information  about  breasts  that  is  never  seen  before.  Yet,  at  its  early  stage  of  development, 
breast  CT  technology  need  to  be  advanced. 

What  this  project  is  addressing  is  the  post-acquisition  image  processing  techniques  to 
improve  the  reconstructed  breast  CT  image  quality.  Firstly,  the  raw  projection  images  of 
breast  CT  are  acquired  without  using  an  anti-scatter  grid.  The  residual  scattered  radiation 
in  the  projections  has  to  be  reduced  via  some  post-acquisition  techniques.  Secondly,  the 
breast  CT  data  are  obtained  with  the  same  dose  as  the  standard  two-view  mammography. 
When  this  level  of  dose  is  equally  split  into  around  500  projections,  each  projection 
contains  considerable  noise.  So  does  the  reconstructed  volume.  Therefore,  it  is  desirable 
to  have  some  image  processing  tools  for  noise  removal. 

The  proposed  project  is  a  collaborative  effort  between  our  group  at  Duke  and  Dr. 
Boone’s  group  in  University  of  California  Davis.  Based  on  the  raw  data  provided  by 
them,  we  will  develop  the  techniques  for  image  improvement  via  the  scatter 
compensation  and/or  denoising. 

Report  Body 

Task  1:  Develop  and  test  a  unique  two-dimensional  Bayesian  image  processing 
technique  on  the  projection  data  of  cone-beam  breast  Computed  Tomography 
(breast  CT)  obtained  without  a  grid. 

This  task  has  been  completed  and  the  results  are  incorporated  into  the  papers  listed  under 
the  category  of  reportable  outcomes.  This  task  is  split  into  two  subtasks.  The  first  subtask 
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is  to  develop  an  algorithm  for  scattered  radiation  removal.  The  second  subtask  is  to 
reduce  the  quantum  noise  from  breast  CT  data. 

Firstly,  an  algorithm  is  designed  and  implemented  for  scatter  reduction.  To  account  for 
the  energy- integrating  characteristic  of  the  flat-panel  digital  detector,  Gaussian 
distributions  are  used  to  approximate  the  signals  recorded  at  individual  pixels  of  the 
detector.  For  the  task  of  removing  scattered  radiation,  the  Gaussian  noise  model  is 
proposed.  The  Maximum  Likelihood  Estimator  (MLE)  is  obtained  via  an  Expectation 
Maximization  algorithm  in  an  iterative  manner.  As  the  image  is  processed  through  many 
iteration  steps,  the  high  frequency  noise  in  the  image  will  be  also  amplified.  In  order  to 
suppress  this  side  effect,  the  Maximum  A  Posteriori  (MAP)  estimator  is  obtained  by 
combining  the  Gaussian  noise  model  with  a  Gibbs  prior  via  Bayes  rule.  It  is  calculated  by 
the  procedure  proposed  by  Hebert  and  Leahy11. 

Figure  1  shows  the  comparison  between  the  original  image  and  MAP  estimate  of  the 
scatter-free  image.  Table  1  is  the  residual  scatter  fraction  (RSF)  and  contrast  to  noise 
ratio  (CNR)  measurements  on  these  three  images.  It  is  shown  that  with  our  algorithm,  the 
scattered  radiation  on  the  images  acquired  without  a  grid  can  be  reduced  down  to  the 
level  achieved  by  using  a  grid.  Meanwhile,  the  CNR  of  the  processed  image  is  twice  that 
of  the  image  acquired  with  a  grid. 


(a)  (b)  (c) 


Figure  1 :  Radiographs  of  an  anthropomorphic  breast  phantom  acquired  on  Siemens  prototype  FFDM 
system  (a)  with  and  (b)  without  an  anti-scatter  grid;  ( b )  the  MAP  algorithm  processed  image  based  on  (b). 
The  white  disks  are  the  beam  stop  (made  of  lead)  array  for  scatter  radiation  measurement. 


Table  1:  Corresponding  residual  scatter  fraction  (RSF)  and  contrast  to  noise  ratio  (CNR)  results  for  the 
three  images  shown  in  Figure  1. 


With  grid 

Without  grid 

Without  grid;  scatter  reduction 

RSF 

11% 

45% 

10% 

CNR 

7.04 

6.99 

15.29 

Secondly,  a  partial-diffusion-equation  based  denoising  technique  was  developed  for  noise 
removal  in  breast  CT.  It  is  observed  that  the  line  integral  images  (converted  from  the  raw 
images  by  logarithmic  operation)  have  the  following  property:  line  integrals  through  the 
center  of  the  breast  have  much  higher  variance  than  those  through  the  edge  of  the  breast. 
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We  derived  a  theoretical  formula  between  the  line  integral  variance  and  the  number  of 

photons  hitting  the  detector  for  a  specific  pixel  region:  var(/ ..)  =  — .  The  details  of  the 

A, 

derivation  can  be  found  in  the  paper  corresponding  to  reportable  outcome  #1.  Based  on 
this  formalism  we  propose  a  PDEtomo  (abbreviation  for:  PDE  for 
tomography/tomosynthesis)  algorithm  for  breast  CT  data. 


Task  2:  Reconstruct  the  three-dimensional  breast  image  based  on  the  processed 
projection  data  from  Task  1. 

This  task  has  already  been  completed  and  some  results  are  shown  in  the  reportable 
outcome  #1.  A  Feldkamp-type  filtered  back  projection  (FBP)  algorithm12  was  custom- 
written  and  used  for  the  cone-beam  reconstruction  of  the  breast  CT  data. 

For  the  results  shown  in  Figure  2,  the  PDEtomo  technique  was  applied  to  the  line  integral 
images  converted  from  raw  projections.  The  processed  projection  images  were  then  fed 
into  the  FBP  core  for  reconstruction.  The  reconstructed  breast  CT  volume  provides 
unique  anatomic  information  of  the  breast  that  is  not  available  before.  In  addition,  the 
PDEtomo  technique  is  very  effective  in  removing  the  noise  while  maintaining  the  details. 


Original 


PDEtomo 


Figure  2:  Reconstructed  coronal  sections  of  a  breast  of  a  human  subject.  The  section  thickness  is  0.5mm. 
The  top  row  is  derived  with  the  original  dataset.  The  bottom  row  is  derived  with  the  PDEtomo  processed 
dataset.  It  is  manifest  that  PDElomo  processed  volume  has  remarkable  less  noise  than  the  original  volume. 


Task  3:  Apply  the  algorithm  in  Task  1  to  the  two-dimensional  slices  of  the 
reconstructed  three-dimensional  breast  image  from  the  unprocessed  projection  data. 

This  task  has  been  completed.  A  two-dimensional  PDE  denoising  technique  (denoted  by 
PDE2d)  was  developed  and  coded. 
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The  PDE2d  algorithm  was  applied  after  reconstruction  on  each  of  the  reconstructed  slices 
to  remove  the  noise.  By  contrast,  the  PDEtomo  technique  presented  in  Task  1  is 
intrinsically  a  three-dimensional  technique,  which  is  applied  to  the  projection  images 
before  tomographic  reconstruction.  The  results  of  these  two  techniques  are  compared 
with  each  other.  One  such  example  is  shown  in  Figure  3.  Figure  3(a)  is  the  same  slice  as 
shown  in  #2  in  the  bottom  row  of  Figure  2.  That  is,  the  volume  was  processed  by  the 
PDEtomo  algorithm.  Figure  3(b)  is  the  one  processed  by  PDE2d  algorithm.  The  parameters 
of  the  techniques  are  adjusted  such  that  they  have  about  the  same  level  of  noise  within  the 
fatty  tissue  region.  However,  Figure  3(a)  has  much  better  anatomical  depiction  than 
Figure  3(b).  The  likely  explanations  for  this  are  as  follows.  Because  the  existence  of  the 
noise  (both  quantum  noise  and  electronic  noise)  in  the  projection  images,  some 
anatomical  details  are  lost  during  tomographic  reconstruction.  By  removing  the  noise  first, 
as  the  PDEtomo  technique  does,  some  anatomical  details  can  be  recovered.  By  contrast, 
since  PDE2d  technique  is  used  after  FBP  reconstruction,  it  cannot  restore  the  lost 
anatomical  detail  information.  Thus,  we  concluded  that  the  denoising  strategy  in  this 
specific  task  is  inferior  to  the  one  in  Task  1. 

(a)  (b) 


Figure  3:  Processed  and  reconstructed  coronal  sections  of  the  same  breast  as  in  Figure  2,  using  the 
denoising  strategies  in  (a)  Taskl  and  (b)  Task3,  respectively.  The  two  slices  have  approximately  the  same 
level  of  noise  removal.  The  one  on  the  left  shows  more  anatomical  details  than  the  one  on  the  right. 

Task  4:  Develop  and  test  three-dimensional  Bayesian  image-processing  technique  on 
the  reconstructed  image  based  on  the  unprocessed  projection  data  acquired  without 
a  grid. 

This  task  has  almost  been  completed.  A  three-dimensional  PDE  technique  (PDE3d)  has 
been  developed  and  is  currently  under  further  evaluation.  The  PDE3d  technique  differs 
from  PDE2d  in  that  it  has  a  three-dimensional  neighborhood  system. 

In  addition,  we  are  exploring  a  PDE3d  counterpart  to  be  applied  before  tomographic 
reconstruction. 

This  task  should  be  completed  within  half  a  month. 
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Task  5:  Develop  a  Computer  Aided  Diagnosis  tool  for  detecting  breast  mass  lesions 
based  on  the  projection  data. 

This  task  is  a  work  in  progress.  Some  preliminary  work  has  been  done. 

The  20  human  volunteer  datasets  were  requested  from  the  collaborating  group  in 
University  of  California  Davis.  The  breast  volumes  were  reconstructed  for  each  case. 
Then  the  afore-mentioned  PDEtomo  technique  was  applied  to  the  breast  CT  datasets  as 
described  in  Task2. 

A  computer-based  virtual  x-ray  CT  imager  is  developed.  It  uses  the  same  parameters  as 
the  physical  breast  CT  system.  A  routine  for  simulation  of  masses  with  random  shape  is 
also  written.  Putting  the  simulated  mass  into  the  virtual  imager,  the  projection  images  can 
be  obtained.  The  simulated  projection  images  with  known  mass  are  then  obtained  by 
summing  up  the  projection  images  from  the  virtual  imager  and  the  clinical  projection 
images. 

The  CAD  tool  based  on  signal  known  exactly  (SKE)  scenario  is  under  development. 

Task  6:  Test  and  compare  the  performances  of  the  CAD  developed  in  Task  5 
applied  to  processed  projection  data  from  Task  1  with  the  CAD  performance  on  the 
projection  data  without  Bayesian  processing. 

Once  the  CAD  in  task  5  is  completely  developed,  this  task  will  be  carried  out  to  compare 
the  CAD  performances  on  original  and  processed  datasets. 


Key  Research  Accomplishments 

•  Proposed  and  developed  the  Gaussian  noise  model; 

•  Proposed  and  developed  a  PDEtomo  algorithm  for  volume  denoising  in  breast  CT ; 

•  Clinical  breast  CT  data  have  been  reconstructed; 

•  Application  of  the  PDEtomo  algorithm  on  simulation  data  and  clinical  data  has 
shown  improvement  of  image  quality  using  CNR  and  resolution  as  figures  of 
merit. 

•  Developed  PDE2d  and  PDE3d  techniques; 

•  A  virtual  breast  CT  imager  was  developed. 

•  A  routine  for  3D  mass  simulation  with  random  shape  was  written. 

Reportable  Outcomes 

Two  publications  were  resulted  from  this  project  in  the  second  fiscal  year,  including  a 
peer-reviewed  paper  and  a  conference  proceeding.  They  are  attached  to  the  end  of  the 
report  as  appendices.  Both  the  name  of  PI  (Jessie  Q.  Xia)  and  the  name  of  the  advisor 
(Joseph  Y.  Lo)  are  bolded. 
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1.  Jessie  Q.  Xia,  Joseph  Y.  Lo,  Kai  Yang,  Carey  E.  Floyd  Jr.,  John  M.  Boone,  Dedicated 
Breast  CT:  Volume  Image  Denoising  via  a  Partial  Diffusion  Equation  Based  Technique. 
Medical  Physics  (submitted). 

2.  Jessie  Q.  Xia,  Georgia  D.  Tourassi,  Joseph  Y.  Lo,  Carey  E.  Floyd  Jr.,  On  the 
Development  of  a  Gaussian  Noise  Model  for  Scatter  Compensation.  Proceedings  of  S PIE 
Medical  Imaging  2007  (in  press). 


Conclusions 

Dedicated  breast  CT  imaging  is  a  promising  technique  for  breast  cancer  imaging.  Since  it 
can  totally  remove  the  overlapping  of  tissues,  it  will  be  even  more  beneficial  for  women 
with  dense  breasts. 

Still,  there  is  a  lot  to  be  done  for  advancing  the  breast  CT  technology.  One  direction  of 
development  is  to  improve  the  image  quality  via  some  post-acquisition  processing 
techniques,  which  is  the  goal  of  this  project.  Specifically,  two  subtasks  are  considered:  1) 
to  remove  scattered  radiation,  and  2)  to  remove  noise.  So  far,  some  techniques  have  been 
designed  and  developed,  and  results  show  that  they  are  effective  for  the  specific  tasks. 

The  project  has  progressed  to  the  stage  where  a  computer  aided  diagnosis  tool  is  under 
development  for  the  assessment  of  mass  detectability  based  on  the  original  breast  CT 
volumes  and  those  with  image  processing. 
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Abstract: 

Dedicated  breast  CT  imaging  possesses  the  potential  for  improved  lesion  detection  over 
conventional  mammograms,  especially  for  women  with  dense  breasts.  The  breast  CT 
20  images  are  acquired  with  a  glandular  dose  comparable  to  that  of  standard  two-view 
mammography  for  a  single  breast.  Due  to  dose  constraints,  the  reconstructed  volume  has 
a  non-negligible  quantum  noise  when  thin  section  CT  slices  are  visualized.  It  is  thus 
desirable  to  reduce  noise  in  the  reconstructed  breast  volume  without  loss  of  spatial 
resolution.  In  this  study,  a  partial-diffusion  equation  (PDE)  based  denoising  technique 
25  specifically  for  breast  CT  was  developed  and  applied  on  the  projection  data  prior  to 
reconstruction.  Simulation  results  show  that  the  PDE  technique  outperforms  Wiener 
denoising.  At  the  photon  fluence  level  of  2.5e4  for  each  projection,  the  noise  of  PDE 
denoised  image  was  39.3%  of  Wiener  denoised  images,  while  the  resolution  was  higher. 


The  PDE  technique  increases  its  performance  advantage  relative  to  Wiener  techniques 


30  when  the  photon  fluence  is  reduced.  For  subjective  evaluation,  the  PDE  technique  was 
applied  to  two  human  subject  breast  datasets  acquired  on  a  prototype  breast  CT  system. 
The  denoised  images  had  appealing  visual  characteristics  with  much  lower  noise  levels 
and  improved  tissue  textures  while  maintaining  sharpness  of  the  original  reconstructed 
volume. 

35 
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Introduction 


40  The  most  common  cancer  type  that  affects  women  globally  other  than  skin  cancer  is 
breast  cancer  Moreover  breast  cancer  is  a  leading  cause  of  cancer- related  women 
mortality,  secondary  only  to  lung  cancer.  It  is  estimated  that  the  disease  will  kill  about 
40,460  US  women  in  2007  '.  Although  mammography  is  the  standard  clinical  screening 
technique  2  3  for  breast  imaging,  superimposition  of  normal  anatomical  structures  may 
45  potentially  obscure  a  breast  lesion.  The  situation  gets  even  worse  for  women  with  dense 
breasts  4,  which  have  more  anatomical  noise  in  the  projection  image.  Researchers  are 
developing  alternative  x-ray  breast  imaging  techniques  that  may  overcome  the  limitations 
of  mammography,  including  three-dimensional  imaging  techniques  such  as  breast 
tomosynthesis  5  6  and  dedicated  breast  CT  7  '2. 

50  Not  long  after  the  CT  technique  was  invented  in  1972,  a  group  of  researchers  studied 
breast  CT  imaging  1 3.  They  applied  the  whole-torso-scanning  mode  and  found  that  a  high 
patient  dose  was  needed  to  achieve  adequate  image  quality.  With  the  advent  of  high- 
resolution  flat-panel  detectors  at  the  end  of  the  1990s,  the  concept  of  breast  CT  once 
again  came  into  researchers’  horizon.  In  particular,  a  2001  paper  7  showed  that  dedicated 
55  breast  CT  could  achieve  quality  breast  images  with  dose  levels  comparable  with  two- 
view  mammography  for  the  same  breast. 

Preliminary  human  subject  data  acquired  on  our  first  prototype  breast  CT  system  14 
provide  exciting  new  information  of  the  breast  that  was  not  available  in  the  past. 
However,  because  the  relatively  low  total  dose  must  be  split  among  a  large  number  of 
60  projection  views  (around  500),  reconstructed  breast  CT  thin  sections  contain  considerable 
quantum  noise.  Thus  it  is  desirable  to  reduce  noise  levels  to  improve  the  conspicuity  of 


breast  lesions.  At  the  same  time,  it  is  desirable  to  retain  spatial  resolution.  Alternatively, 
by  applying  denoising  techniques,  the  gain  of  improving  CT  image  quality  can  be 
exchanged  for  dose  reduction  while  maintaining  the  image  quality. 

65  For  low  dose  CT,  some  general-purpose  sinogram  smoothing  techniques  based  on  either 
penalized  likelihood  15  or  penalized  weighted  least  squares  16  were  developed.  These 
techniques  can  be  potentially  applied  on  dedicated  breast  CT  datasets.  Zhong  et  al  17 
developed  a  wavelet-based  technique  and  applied  it  on  phantom  breast  CT  data.  Their 
results  showed  that  with  denoising,  dose  could  be  potentially  reduced  by  up  to  60%. 

70  The  Partial  Diffusion  Equation  (PDE)  based  technique  18' 19  is  another  state-of-the-art 
denoising  method  which  is  effective  not  only  in  removing  noise  but  also  in  preserving 
detail.  Although  computationally  intensive,  this  iterative  method  can  provide  more 
freedom  in  choosing  the  desired  denoising  effect.  In  this  study,  we  describe  a  PDE  based 
denoising  technique  specifically  for  breast  CT  and  evaluate  it  both  on  simulated  and 

75  empirically  collected  human  subject  datasets. 

The  image  quality  of  PDE  denoised  images  was  compared  against  that  of  standard 
Weiner  filtering  techniques.  Quantitative  comparisons  were  made  using  simulated  data  at 
various  exposure  levels,  while  qualitative  comparisons  were  made  using  dedicated  breast 
CT  scan  data  from  two  human  subjects. 


80 


Materials  and  Methods 


85  Dedicated  Breast  CT  System  and  Human  Subject  Datasets 

As  is  illustrated  in  Figure  1,  dedicated  breast  CT  systems  are  typically  designed  as 
follows:  a  patient  lies  supine  on  a  lead-shielded  table  with  one  breast  hanging  through  a 
hole  on  the  table  in  the  pendant  geometry.  The  x-ray  tube  and  the  flat  panel  detector 
rotate  in  the  horizontal  plance  underneath  the  table.  This  setup  is  different  from  a 
90  conventional  CT  system,  where  the  x-ray  tube  and  detector  rotate  around  the  torso  of  a 
patient  (axial  scanning).  Since  only  the  breast  to  be  imaged  is  exposed  to  the  x-ray  beam, 
the  dose  to  the  patient  can  be  greatly  reduced.  A  pilot  study  7  showed  that  this  type  of 
dedicated  breast  CT  system  is  able  to  achieve  a  satisfactory  image  quality  with  dose 
levels  comparable  to  standard  two- view  mammography  for  the  same  breast. 

95  Using  the  above  system  design,  a  custom-designed  breast  CT  system  was  fabricated  at 
the  University  of  California  Davis  Medical  Center,  and  is  currently  accruing  patient 
images.  The  x-ray  tube  has  a  Comet  beryllium- windowed,  water-cooled  tungsten  anode 
and  a  nominal  focal  spot  with  the  size  of  0.4  mm  x  0.4  mm.  A  Pantak  high  frequency  x- 
ray  generator  drives  the  tube  with  the  voltage  ranging  from  10  kV  to  110  kV.  The  Csl- 
100  based  flat  panel  detector  (Varian,  PaxScan  4030CB)  has  a  field  of  view  of  40  cm  x  30  cm. 
Using  30  frames  per  second  and  2x2  pixel  binning  mode,  the  detector  generates  the 
images  each  with  matrix  size  of  1024  x  768  and  pixel  dimension  of  0.388  mm  x  0.388 
mm.  A  Kollmorgen  servo  motor  was  employed  to  drive  the  rotation  of  the  tube-detector 
gantry  as  well  as  encode  the  angular  information.  The  source-to-isocenter  distance  is  46.9 
105  cm  and  the  source-to- detector  distance  is  88.4  cm. 


For  the  two  human  subject  datasets  presented  in  this  paper,  the  projection  images  were 
acquired  under  80  kVp  using  a  circular  orbit.  The  mAs  values  were  chosen  for  each 
subject  such  that  the  mean  glandular  dose  using  breast  CT  was  equal  to  two-view 
mammography.  Each  subject  is  scanned  within  17  seconds  to  get  a  total  of  500  projection 
110  images  that  span  slightly  over  360  degrees.  After  dead  pixel  and  flat  field  corrections, 
each  dataset  is  ready  for  tomographic  reconstruction. 

Simulated  Breast  CT  Datasets 

In  this  study,  some  simulated  breast  CT  datasets  were  also  generated  to  aid  the  analysis. 
The  computer-generated  breast  is  a  hemisphere  with  radius  of  7cm.  It  has  homogeneous 
115  breast  tissue  with  a  uniform  linear  attenuation  coefficient  of  0.28  cm  1  and  is  surrounded 
by  1mm  thick  skin  20  with  linear  attenuation  coefficient  of  0.5  cm1.  A  spherical  lesion 
with  10%  contrast  relative  to  the  breast  tissue  was  embedded  at  the  center  of  the  breast. 
This  simulated  breast  was  scanned  virtually  by  a  monochromatic  x-ray  cone  beam  with 
infinitely  small  focal  spot  and  ideal  flat-panel  detector  with  100%  detective  quantum 
120  efficiency.  The  geometric  parameters  are  the  same  as  the  physical  breast  CT  scanner 
described  in  the  previous  subsection. 

For  each  2D  projection  image,  an  analytical  line  integral  image  was  first  obtained  based 
on  the  aforementioned  virtual  dedicated  breast  CT  scanning.  A  noisy  raw  image  was 
generated  according  to  the  measurement  model  21 .  The  model  takes  into  account  both 
125  photon  quantum  noise  and  electronic  read-out  noise.  It  has  the  following  form: 

Yj  =  Gi  E  ■  Poisson(I0e~l‘ )  +  Gaussian(0,a2) ,  (1) 

where  G,  is  the  gain  factor  of  the  imaging  system,  E  is  the  mean  energy  level  of  the 

polychromatic  x-ray  beam,  and  the  Gaussian  term  is  for  the  electronic  noise.  In  our 


simulation  we  chose  GpO.0035  /kcV,  £'=40  keV,  and  cr=10.  The  values  of  G;  and  cr 
130  were  referred  to  those  used  in  21 . 

Tomographic  Reconstruction 

Before  reconstruction,  the  raw  projection  images  undergo  a  preprocessing  step.  The 
simplest  preprocessing  step  is  to  convert  the  raw  projection  image  into  the  line  integral 
via  the  logarithm  operation.  On  each  raw  projection  image,  a  region  of  interest  (ROI)  is 
135  identified  which  is  outside  the  breast  silhouette,  and  the  pixel  value  without  attenuation 
70  is  approximated  by  the  mean  pixel  value  with  the  ROI.  Then  the  line  integral  image  is 
obtained  by:  /;/  =  log(/0  //i;),  where  ItJ  is  the  pixel  value  at  (i,  j)  position. 

Since  the  datasets  have  high  angular  sampling  rate,  the  computationally  efficient  filtered 
backprojection  (FBP)  22  algorithm  was  chosen  for  reconstruction.  The  Feldkamp  type 
140  FBP  for  cone-beam  geometry  was  custom  written  and  a  Shepp-Logan  filter  was  used. 

PDE  Denoising  Technique 

The  PDE  denoising  strategy  was  inspired  by  the  heat  equation,  which  characterizes  the 
spatial  temperature  distribution  as  a  function  of  time.  When  the  time  period  is  very  long, 
the  temperatures  tend  to  be  uniform  in  space,  i.e.,  the  temperature  image  is  smoothed. 
145  Therefore,  if  an  image  needs  to  be  denoised,  it  can  be  treated  as  a  temperature  image  and 
fed  into  the  heat  equation  for  its  smoothed  versions  as  time  goes  by.  The  simplest  case 
will  be  the  heat  propagation  in  a  homogenous  media,  or  equivalently,  smoothing  equally 
all  over  the  space.  The  heat  equation  in  this  case  will  be: 

—  =  kV2/  =  V-(kV/),  (2) 

dt 

150  where  k  is  a  scalar  diffusivity  constant  of  the  homogeneous  media,  V/  is  the  gradient  of 
the  image  I  and  V2/  =  V  •  (V/)  is  the  Laplace  operation  on  I  over  the  spatial  variables  23. 


This  simplest  case  will  smooth  out  both  the  noise  and  the  detail  such  as  the  edges  in  the 
image.  However,  a  nonlinear  PDE  can  reduce  noise  while  preserving  spatial  resolution  in 
the  image: 

155  ^-V-QK|V/|)W),  (3) 

at 

where  p(.)  is  called  the  diffusivity  function,  a  function  of  the  gradient  norm  I VII.  Equation 
(2)  is  a  special  case  of  Equation  (3)  where  p(.)  is  chosen  to  be  a  constant  value  k.  The  p(.) 
named  diffusivity  is  a  function  to  regulate  the  local  smoothness.  The  choice  of  p(.)  is 
critical  to  PDE-based  denoising  of  images.  In  this  work,  we  chose  a  diffusivity  function 
160  proposed  by  Perona  and  Malik  18: 

_d^_ 

p(d)  =  e~s 2  (4) 

where  <5  is  a  user- specified  parameter.  When  the  image  gradient  norm  is  very  large  at  a 
location  region,  the  diffusivity  will  be  very  small,  and  thus  the  local  image  values  will  be 
preserved  within  a  small  time  period  whereas  another  more  uniform  region  will  be 
165  smoothed  out  at  the  same  time.  The  parameter  8  acts  like  a  cut-off  value;  image  regions 
with  gradient  norm  below  5  will  have  more  noise  removed  while  regions  with  a  higher 
gradient  norm  will  stay  sharp. 

The  diffusion  equation  can  be  discretized  by  the  finite  difference  approach  using  the  first- 
order  neighborhood  system.  Each  pixel  has  four  neighbors:  the  north,  south,  west  and  east 
170  neighbor  pixels.  Assuming  Ax  =  Ay  =  1  in  the  two-dimensional  case,  the  discretized 
version  of  Equation  (3)  is 


i.'*' 

"  .  "  =  IPn  ■  W7  +  Ps  ■  VT  +  Pw-VwI+Pe-  V£/]7 


(5) 


where  A t  is  discretized  time  step,  [pN,  ps,  pw,  pE  ]  are  the  values  of  the  diffusivity  function 
at  the  neighboring  pixels  of  location/^ ,  and  V*I  is  a  notation  for  the  difference  between 

175  the  pixel  value  of  one  neighbor  and  itself.  In  this  work.  At  =  0.1  was  used. 

The  simulated  breast  CT  data  revealed  a  pattern  of  noise  on  its  line  integral  images.  That 
is,  toward  the  chest  wall  of  the  breast  the  noise  is  larger.  When  the  photon  fluence  is 
reduced,  the  phenomenon  becomes  even  more  obvious.  A  line  profile  is  shown  in  Figure 
2  to  help  illustrate  this  point.  It  can  be  explained  theoretically.  Again,  a  simplifying 
180  assumption  of  monochromatic  beam  is  used.  If 
/,.  ~  Poissoni A.  ) 

y  y 

,  ,  ,  then  the  variance  of  the  line  integral/,,  can  be  approximated 

/,/='ogTL  =  'og/0-logA  y 

y 

by  the  delta  method  24  using  the  second-order  Taylor  expansion: 

var (ly)  =  var(log/.  )  =  var (/..)  *  [/'(log/,7)]2  a  (6) 

'ti 

This  formalism  can  be  integrated  into  the  PDE  denoising  technique  by  adapting  the 
185  parameter  5  in  the  diffusivity  function  spatially  as: 


The  resultant  PDE  denoising  technique  is  denoted  PDE  for  Tomography/Tomosynthesis 
(PDEtomo). 

The  PDE  denoising  method  can  be  applied  to  the  CT  data  at  different  points  in  the 
190  reconstruction  process.  When  a  Filtered  backprojection  (FBP)  algorithm  is  used  for 
reconstruction,  there  are  four  potential  locations  where  the  denoising  techniques  can 
applied,  which  are  illustrated  in  Figure  3.  In  this  study,  the  denoising  techniques  were 


applied  at  step  2.  That  is,  the  line  integral  images  are  processed  and  these  denoised  data 
were  used  in  the  FBP  algorithm. 

195  The  datasets  were  also  processed  by  Wiener  adaptive  filtering  technique  commercially 
implemented  by  MATLAB  (Mathworks  Inc,  Natick,  Massachusetts),  for  comparison. 

Image  Evaluation 

Contrast  of  the  lesion,  noise  level,  contrast-to-noise  ratio  (CNR)  and  resolution  25  are 
metrics  for  the  quantitative  evaluation  of  the  denoising  technique.  These  metrics  were 
200  calculated  based  on  the  reconstructed  slices  of  the  simulated  breast.  The  contrast  of  the 
lesion  is  defined  as  the  relative  difference  between  the  average  pixel  values  within  the 
lesion  and  those  outside  the  lesion.  The  noise  level  is  characterized  by  the  percent  noise, 
that  is,  the  standard  deviation  of  the  pixel  values  within  a  uniform  ROI  relative  to  the 
mean  value.  This  is  also  called  the  coefficient  of  variation  (COV).  The  CNR  measure  is 
205  the  ratio  of  contrast  of  the  lesion  to  the  percentage  noise,  which  is  equivalent  to  the  signal 
to  noise  ratio  (SNR). 

In  order  to  measure  spatial  resolution,  a  high  intensity  sphere  was  simulated  within  the 
breast,  projected,  and  target  reconstructed  (with  an  in-plane  pixel  dimension  of  0.2mm) 
around  the  high-intensity  sphere.  Then  the  edges  of  the  circular  disc  within  a 
210  reconstructed  slice  are  averaged  radically.  For  each  specific  denoising  technique,  a 
Gaussian  function  is  fitted  to  the  edge,  and  the  full  width  at  half  maximum  (FWHM)  of 
that  function  was  measured  as  the  parameter  to  evaluate  the  spatial  resolution. 


Results 


Simulation  Results 

215  In  Figure  4,  the  three  columns  correspond  to  the  reconstructed  slices  of  the  simulated 
breast  CT  data  without  denoising,  with  Wiener  denoising  and  with  the  PDEtomo 
denoising,  respectively.  The  four  rows  show  the  results  using  varying  photon  fluences 
using  I0=2.5e3,  5e3,  7.5e3  and  2.5e4  correspondingly  from  the  top  to  the  bottom.  The 
noise  level  is  appreciably  lower  for  each  successive  row.  The  reconstructed  slice 
220  thickness  was  0.5  mm  and  within-plane  pixel  dimension  was  0.8  mm.  These  slices 
contain  a  lesion  at  the  center,  which  is  only  barely  visible  for  the  highest  fluence  level 
(bottom  panel).  With  Wiener  processing,  the  lesion  is  visible  for  the  highest  two  fluence 
levels.  With  PDEtomo  processing,  the  lesion  is  visible  for  the  highest  three  fluence  levels. 
To  demonstrate  the  effect  of  denoising  on  the  uniformity  of  the  noise  level  across  the 
225  image,  another  reconstructed  coronal  view  slice  containing  only  homogeneous  breast 
tissue  is  shown  in  Figure  5(a).  For  the  horizontal  center  line,  the  standard  deviation  of  a 
local  7x7  ROI  centered  at  each  pixel  on  that  horizontal  line  is  plotted  in  Figure  5(b)  to  (d) 
for  photon  fluence  of  2.5e3,  5e3  and  2.5e4. 

In  all  three  cases,  the  original  image,  as  shown  on  the  left  of  Figure  5(a),  presents  a  non- 
230  uniform  level  of  noise  such  that  the  interior  section  has  higher  standard  deviation  values, 
that  is,  higher  noise  levels  than  the  periphery.  Wiener  processing  reduces  the  overall 
noise  level,  but  the  non-uniform  trend  remains  in  Figure  5(b)  and  6(c)  corresponding  to 
lower  photon  fluences.  The  PDEtomo  processed  images  reduce  the  noise  and  also  the 
non-uniformity  of  the  plot.  For  the  highest  fluence  images  in  Figure  5(d),  the  effects  of 
235  Wiener  and  PDEtomo  denoising  are  very  similar. 


The  CNR  results  based  on  the  images  depicted  in  Figure  4  are  given  in  Figure  6.  At  the 
highest  fluence  of  2.5e4,  CNRs  of  the  Wiener  and  PDEtomo  processed  images  are  3.59 
and  4.08  times  that  of  the  original  one,  respectively.  The  paired  two-tailed  student  t  test 
of  CNRs  of  the  Wiener  and  PDEtomo  results  gives  the  p-value  of  0.0036,  indicating  that 
240  the  CNRs  of  PDEtomo  were  statistically  higher  than  those  of  Wiener.  Similar  trends  were 
observed  at  the  lower  fluence  levels. 

The  resolution  results  are  given  in  Figure  7.  The  first  image  in  the  top  panel  shows  the 
full  view  reconstructed  slice  with  a  high-intensity  sphere  embedded  into  the  center  of  the 
breast  tissue.  The  ROI  for  targeted  reconstruction  is  marked  on  the  slice.  The  original 
245  reconstructed  slice  as  well  as  the  wiener  denoised  one  are  to  the  right  of  the  first  row  of 
Figure  7(a). 

In  the  second  row,  the  PDEtomo  denoised  result  is  shown  in  the  middle.  To  its  left,  a 
Gaussian  kernel  is  directly  applied  to  the  original  target-reconstructed  slice  and  its 
resultant  blurred  slice  obtains  the  same  noise  level  as  the  PDEtomo  result.  To  its  right,  a 
250  corresponding  Gaussian  blurred  one  obtains  the  same  resolution  level.  As  a  comparison, 
the  PDEtomo  achieves  a  low  noise  level  and  the  high  resolution  simultaneously.  In 
Figure  7(b),  the  noise  and  resolution  values  are  obtained  and  plotted  for  both  the 
PDEtomo  and  Wiener  techniques.  As  shown  in  the  noise-resolution  plot,  PDEtomo 
processing  outperforms  Wiener  processing  because  it  has  lower  noise  and  higher 
255  resolution.  The  noise  levels  of  the  Wiener  and  PDEtomo  techniques  are  12.0%  and  4.7%, 
respectively.  The  noise  level  of  Wiener  is  2.54  times  that  of  the  PDEtomo.  At  the  same 
time,  the  resolution  measure  was  0.2209  mm  for  PDEtomo  and  0.2867  mm  for  Wiener. 


The  resolution  of  PDEtomo  is  better  than  that  of  Wiener. 
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Human  Subject  Results 


Figure  8  and  9  show  the  coronal  reconstructed  slices  from  two  human  subject  breast  CT 
data.  By  visual  comparison,  PDEtomo  technique  (bottom  row)  reduces  the  noise 
considerably  while  maintaining  the  resolution  of  the  original  reconstruction  (top  row). 
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Discussion 


Dedicated  breast  CT  imaging  is  an  exciting  new  candidate  that  possesses  the  potential  of 
improving  breast  cancer  diagnosis  over  conventional  mammography.  Some  preliminary 
studies  14' 26' 27  show  that  satisfactory  images  can  be  acquired  using  the  same  dose  as 
270  standard  two-view  mammography.  When  this  dose  is  divided  among  the  potentially 
hundreds  of  individual  projection  images  comprising  each  scan,  the  high  level  of 
quantum  noise  in  the  projection  data  will  pass  through  to  the  final  reconstructed  volume. 
This  motivates  the  development  of  denoising  tools  to  effectively  remove  the  noise, 
improve  lesion  conspicuity,  and  maintain  image  resolution. 

275  In  this  work,  a  partial  diffusion  equation  (PDE)  based  denoising  technique  was 
specifically  developed  for  processing  breast  tomography  data.  This  PDEtomo  technique 
takes  into  account  the  noise  distribution  characteristic  in  the  projection  image  after 
converting  to  the  line  integrals  via  the  nonlinear  logarithm  operation. 

Both  the  quantitative  results  in  the  simulation  study  and  the  visual  inspection  in  human 
280  subject  data  study  showed  the  promise  of  this  new  PDEtomo  technique.  In  the  simulation 
study,  it  was  compared  with  Wiener  technique,  an  adaptive  technique  that  is  accepted  as  a 
competitive  denoising  tool.  The  results  show  that  the  new  denoising  technique  can 
achieve  lower  noise  level  in  the  reconstructed  volume  with  higher  resolution  than  Wiener 
technique.  The  low  contrast  lesion  put  in  the  center  of  the  simulated  breast  can  be  better 
285  detected  on  the  PDEtomo  processed  datasets  than  on  the  Wiener  processed  datasets,  as  is 
shown  in  Figure  4.  Not  only  is  the  noise  in  the  reconstructed  slices  filtered  by  PDEtomo 
lower  than  those  by  Wiener  filtering,  but  the  noise  levels  all  over  the  breast  region  tend  to 
be  more  uniform  as  well.  This  is  shown  in  Figure  5.  The  advantage  of  PDEtomo  over 


Wiener  filtering  in  terms  of  decreased  noise  level  and  improved  noise  uniformity  are  both 
290  more  evident  for  lower  dose  cases,  which  indicates  that  the  PDEtomo  technique  holds 
promise  for  processing  datasets  acquired  at  lower  dose  levels. 

Even  though  the  theoretical  description  of  the  noise  variance  in  the  projection  image  due 
to  the  quantum  noise  and  the  logarithm  operation  is  more  approximate  for  the  empirical 
data,  the  PDEtomo  technique  still  provides  good  denoised  images. 

295  There  are  some  limitations  in  this  study.  Firstly,  the  simulated  breast  CT  data  are  based 
on  a  monochromatic  x-ray  beam  with  the  kV  value  set  to  be  approximately  the  same  as 
the  effective  kV  value  of  the  x-ray  beam  used  to  acquire  the  empirical  data.  Secondly,  the 
parameters  in  the  measurement  model  used  for  adding  noise  to  the  simulated  projection 
images  are  all  hypothetical,  given  that  presently  their  empirical  values  are  lacking.  Hence, 
300  the  task  of  calibrating  the  dose  in  the  simulation  study  cannot  be  fulfilled  at  this  stage.  In 
future  work,  considerable  optimization  remains  to  be  performed  to  calibrate  the  PDEtomo 
technique  using  empirical  images  taken  with  physical  phantoms  as  well  as  human 
subjects.  Given  the  robust  trends  shown  in  this  study,  however,  the  PDEtomo  technique 
should  continue  to  match  or  outperform  the  Wiener  technique,  especially  if  dose  is 
305  further  lowered  such  as  to  achieve  a  breast  CT  scan  with  equal  or  less  dose  than  two- view 
conventional  mammograms. 

Due  to  the  very  low  photon  fluence  on  each  projection  view  in  dedicated  breast  CT,  the 
electronic  noise  is  one  of  the  major  sources  of  the  overall  noise,  especially  in  dense  breast 
regions  or  if  the  dose  is  further  reduced.  The  present  version  of  the  PDEtomo  technique 
310  doesn’t  consider  the  effects  of  additive  electronic  noise.  It  will  be  worthwhile  to  explore 
the  possibility  of  taking  the  characteristics  of  this  type  of  noise  into  account  in  the 


denoising  technique  or  to  combine  it  with  a  statistical-modeling  approach  that  explicitly 
treats  the  electronic  noise. 

In  conclusion,  a  partial  diffusion  equation  based  denoising  technique  was  developed 
315  specifically  for  dedicated  breast  CT  data.  By  incorporating  into  the  algorithm  the 
knowledge  of  the  non-uniform  distribution  of  the  noise  in  the  projection  image  after  the 
preprocessing  step,  it  provides  excellent  denoised  data  with  sharp  edges.  The  technique 
shows  the  most  promise  on  datasets  acquired  with  lower  dose. 
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Figure  1  Illustration  of  a  dedicated  breast  CT  system. 

The  x-ray  tube  and  flat-panel  detector  rotate  together  around  the  breast.  And  the 
breast  is  the  only  region  to  be  illuminated. 
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Figure  2  One-dimensional  line  integral  profiles  across  the  breast  on  a  projection 
350  image. 

The  dashed  (blue)  and  continuous  (red)  plots  correspond  to  noise  free  case  and  the 
case  with  10=2. 5e3,  respectively.  The  variance  of  line  integral  is  larger  at  the  center 
of  breast  region,  and  gets  lower  toward  the  skin. 
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Figure  3  Illustration  of  the  possible  locations  of  image  denoising  module  to  be 
applied  with  respect  to  the  reconstruction  module  for  dedicated  breast  CT  data. 

In  this  study,  the  denoising  techniques  are  applied  at  location  2. 
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Figure  4  The  reconstructed  slices  of  original,  Wiener  and  PDEtomo  processed 
datasets  using  the  varying  photon  fluence  levels  from  the  simulation  study. 

The  central  low  contrast  lesion  is  visible  only  at  the  highest  fluence  level  with  the 
original  dataset.  It  is  visible  at  the  highest  two  fluence  levels  with  the  Wiener 
365  processed  dataset,  whereas  it  is  visible  at  the  highest  three  fluence  levels  with  the 
PDEtomo  technique. 
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Figure  5  (a)  The  reconstruction  slices  without  a  lesion.  Along  the  horizontal  lines 
marked  in  (a),  local  noise  level  profiles  are  computed  for  original,  Wiener  and 
PDEtomo  processed  datasets  at  10  of  (b)  2.5e3,  (c)  5e3  and  (d)  2.5e4. 

380  PDEtomo  processed  results  have  lower  and  more  uniform  local  noise  levels  than 
Wiener  results.  The  difference  is  more  manifest  at  lower  fluence  level. 


4.0 

3.5 

3.0 

2.5 

0£ 

Z  2.0 

<J 

1.5 

1.0 

0.5 

0.0 

orig  Wiener  PDE  tomo 


I 


Figure  6  CNR  comparison  between  the  original,  Wiener  denoised  and  PDE  in  a 
simulation  study  with  I0=2.5e4. 

The  PDEtomo  has  consistently  higher  CNRs  than  Wiener  for  the  same  data.  The 
paired  two-tailed  student  t  test  gives  a  p-value  of  0.0036,  which  indicate  the 
difference  between  the  CNRs  of  PDEtomo  and  Wiener  are  statistically  significant. 
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395  Figure  7  (a)  The  target  reconstruction  around  a  high  contrast  object  from  original, 
Wiener  processed  and  PDEtomo  processed  datasets,  and  (b)  noise-resolution  plot 
with  I0=2.5e4. 

The  data  processed  by  PDEtomo  has  a  lower  noise  level  and  higher  resolution  than 
the  one  processed  by  Wiener.  The  two  Gaussian  kernel  blurred  ones  are  to  match 
400  the  noise  level  or  the  resolution  with  the  PDEtomo  processed  data. 


405  Figure  8  Human  Subject  Result  No.l 

Top  row  shows  original  reconstruction,  coronal  slices  from  normal  breast.  Bottom 
row  shows  same  slices  with  PDEtomo  denoising,  resulting  in  remarkably  reduced 
noise  levels. 
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Figure  9  Human  Subject  Result  No.2 

Top  row  shows  original  reconstruction,  coronal  slices  from  normal  breast.  Bottom 
row  shows  same  slices  with  PDEtomo  denoising,  resulting  in  markedly  reduced 
noise  levels. 
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ABSTRACT 

The  underlying  mechanism  in  projection  radiography  as  well  as  in  computed  tomography  (CT)  is 
the  accumulative  attenuation  of  a  pencil  x-ray  beam  along  a  straight  line.  However,  when  a 
portion  of  photons  is  deviated  from  their  original  path  by  scattering,  it  is  not  valid  to  assume  that 
these  photons  are  the  survival  photons  along  the  lines  connecting  the  x-ray  source  and  the 
individual  locations  where  they  are  detected.  Since  these  photons  do  not  carry  the  correct  spatial 
information,  the  final  image  is  contaminated.  Researchers  are  seeking  techniques  to  reduce 
scattering,  and  hence,  improve  image  quality,  by  scatter  compensation.  Previously,  we  presented 
a  post-acquisition  scatter  compensation  technique  based  on  an  underlying  statistical  model.  We 
used  the  Poisson  noise  model,  which  assumed  that  the  signals  in  the  detector  individually 
followed  the  Poisson  process.  Since  most  x-ray  detectors  are  energy  integrating  rather  than 
photon  counting,  the  Poisson  noise  model  can  be  improved  by  taking  this  property  into  account. 
In  this  study,  we  developed  a  Gaussian  noise  model  by  the  matching-of-the-first-two-moments 
method.  The  Maximum  Likelihood  Estimator  of  the  scatter-free  image  was  derived  via  the 
expectation  maximization  (EM)  technique.  The  maximum  a  posteriori  estimate  was  also 
calculated.  The  Gaussian  noise  model  was  preliminarily  evaluated  on  a  full-field  digital 
mammography  system. 

KEYWORDS:  Scatter  Compensation,  Scatter  Reduction,  Gaussian  Noise  Model, 
Expectation  Maximization 


1.  INTRODUCTION 

Scattered  radiation  degrades  medical  images.  A  recent  Monte  Carlo  study  showed  that  scattered 
radiation  causes  the  drop  of  low-frequency  modulation  transfer  function  (MTF),  changes  the 
shape  of  MTF  and  adds  considerable  noise  to  projection  images  1 .  In  computed  tomography  (CT), 
scattered  radiation  leads  to  cupping  artifacts  on  reconstructed  sections.  Therefore,  removal  of 
scattered  radiation  from  projection  images  is  essential  for  improved  image  quality,  particularly 
for  the  latest  advanced  imaging  radiography  techniques  including  dedicated  breast  CT  and  breast 
tomosynthesis,  which  typically  do  not  use  anti-scatter  grids. 

There  are  two  categories  of  scatter  compensation  techniques:  hardware  based  ones  and  numerical 
compensation  methods.  This  study  uses  a  numerical  compensation  method  to  develop  a  statistical 
scatter  reduction  technique.  Previously,  the  exposure  value  of  each  pixel  recorded  by  a  detector 
was  modeled  by  Poisson  distribution  2.  For  flat-panel  detectors,  which  belong  to  the  type  of 
energy  integrating  rather  than  photon  counting,  the  underlying  statistics  is  a  compound  Poisson 
process 3.  It  may  be  well  approximated  by  Gaussian  distribution. 


In  this  study,  we  will  propose  a  Gaussian  noise  model  for  scatter  reduction,  and  derive  a 
maximum  likelihood  expectation  maximization  (MLE  or  MLEM)  algorithm  and  a  maximum  a 
posteriori  (MAP)  algorithm.  We  will  then  apply  the  algorithms  to  radiographs  acquired  on  FFDM 
for  preliminary  evaluation. 

2.  MATERIALS  AND  METHODS 
2.1  Gaussian  Noise  Model 

In  the  chosen  numerical  scatter  compensation  scheme,  the  projection  image  is  the  sum  of  the 
primary  radiation  and  scattered  radiation.  We  have  modeled  scattered  radiation  as  the  two- 
dimensional  convolution  of  primary  radiation  with  a  scatter  kernel,  which  is  displayed  as  a  double 
exponential  function  (Figure  1). 
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Figure  1:  Schematic  of  a  scatter  kernel  with  a  radically  exponential  shape.  It  has  two  parameters:  full 
width  at  half  maximum  and  magnitude). 


By  the  matching-of-the-first-two-moments  method,  we  approximated  the  energy-integrating 
signal  by  using  a  Gaussian  distribution  and  created  a  scatter  compensation  model  called  the 
Gaussian  noise  model. 


d,  \B,ojf  ~  Gaussianfb^Ojf) 

si  IB, cr,,2  ~  Gaussian((B **P)jtop2) 

y.  =  d,  +st  I  B.  a2 .  a  2  ~  Gaussian(bi  +  (B**P)j,oil2  +cri22) 


where  dt,  st,  and  y,  are  pixel  values  at  location  i  corresponding  to  primary,  scattered  and  total 
radiation,  respectively,  B,-  is  the  expectation  of  d,-,  and  a, \  and  on  are  the  variance  of  pixel  values 
related  to  the  primary  radiation  and  scattered  radiation,  respectively. 


Using  the  expectation  maximization  (EM)  algorithm  shown  in  the  appendix,  the  MLE  of  the  ideal 
scatter-free  image  was  derived  with  analytical  form  shown  in  Equation  (2). 

b(rl)  -b™  +"*  -[yk  ~(bk(n>  +(Bin)  **P)k)]  (2) 

wk  =ol  /(cta2,  +o;2) 

MLE  estimate  is  known  to  increase  high  frequency  image  noise.  To  overcome  this,  some 
constraints  can  be  put  on  the  noise  level  within  the  estimated  B,  in  other  words,  prior 
information  about  B  can  be  provided.  Thus,  by  Bayes’s  Rule, 

p(B  I Y)  oc  p(Y  I  B)p(B) ,  (3) 

where  p(B\Y)  is  the  posterior  joint  distribution  of  B,  given  the  measured  pixel  values  Y  =  {y;; 
i=l,...,N},  p(Y\B)  is  equal  to  the  likelihood  of  B,  and  p(B)  is  the  prior  joint  distribution  of  B  = 
{bi;  i=l,...,N}. 

We  assumed  B  is  a  Markov  random  process.  It  therefore  follows  a  Gibbs  distribution: 

p{B)  =  -e-uwl^  (4) 

K 

where  K  is  a  normalizing  factor  which  is  independent  of  B,  U(B)  is  the  energy  function,  and  /3 
is  a  free  parameter  adjusting  the  relative  weight  of  this  prior  on  the  maximum  a  posteriori 
estimator  of  B.  When  /3  approaches  infinity,  the  MAP  of  B  approaches  the  MLE  of  B. 

The  energy  function  is  the  sum  of  the  potential  function,  i.e., 

V (B)  =  ^  VC(B)  >  (5) 

<€C 


where  C  is  the  set  comprised  of  all  cliques  in  the  image.  In  this  study,  the  Gibbs  prior  is 
defined  over  a  second-order  neighborhood  system  (for  each  pixel,  its  north,  south,  east,  and  west 
neighboring  pixels  plus  its  four  diagonal  neighboring  pixels),  with  each  clique  comprised  of  two 
neighboring  pixels.  There  are  many  forms  of  the  potential  function  VC(B).  We  chose  one 
that  is  adaptive  to  discontinuity  4: 


Vc({b,:bJ})  = 


S  e2  +  (bi-bJ)2 


(6) 


where  i  and  j  are  the  neighboring  pixels  within  the  clique  i~j  and  b,  and  bt  represent  their 
respective  intensities.  5C  is  an  adjustable  parameter  to  regulate  the  cut-off  frequency  of  the 
noise  in  the  image. 

The  MAP  estimate  of  {b;}  was  calculated  through  the  two-step  maximization  procedure  proposed 
by  Hebert  and  Leahy  5. 

2.2  Test  Images 

Images  were  acquired  with  a  Siemens  prototype  digital  mammography  system  (Mammomat 
Novation  DR;  Siemens,  Erlangen,  Germany)  with  70  pm  isotropic  resolution.  Uniform  breast 
phantoms  (CIRS,  Inc.,  Norfolk,  VA)  were  imaged  (28kVp),  with  a  Mo/Mo  target/filter 


combination.  The  phantoms  were  designed  to  be  radiographically  equivalent  to  a  4-cm-thick 
compressed  breast  with  50%  glandular  tissue  density.  A  built-in  square-shaped  dent  in  the  center 
of  the  phantom  mimicked  a  high-contrast  lesion  in  the  digital  mammography  images.  All  images 
were  acquired  without  an  anti-scatter  grid.  For  the  purpose  of  scatter  measurement,  all  images 
were  repeated  with  an  array  of  beam  stops  (lead  discs  3  mm  in  diameter)  superimposed  on  the 
breast  phantoms.  Because  lead  discs  absorb  all  the  primary  radiation,  only  scatter  radiation  can  be 
detected  behind  them  (Figure  2). 


Figure  2:  Radiograph  of  the  tissue  equivalent  slabs.  The  arrays  of  black  disks  are  the  shadow  of  beam 
stops.  The  CNR  values  are  obtained  based  on  the  bright  square  region  of  interest. 


Images  were  then  be  fed  into  the  algorithms  for  processing.  The  results  were  then  evaluated 
through  various  metrics  described  in  the  following  subsection. 

2.3  Image  Evaluation  Metrics 

Three  algorithms  were  employed  to  estimate  the  expected  amount  of  primary  radiation.  Its 
effect  was  measured  by  the  residual  scatter  fraction  (RSF).  At  the  same  time,  we  anticipated 
that  the  contrast-to-noise  ratio  (CNR)  would  be  constrained  or  even  improved  after  image 
processing.  We  implemented  a  metric  to  estimate  post-contrast  CNR.  Finally,  we  monitored 
with  a  test  bar  the  effect  of  the  algorithms  on  the  spatial  resolution  of  the  images. 

2.3.1  Residual  Scatter  Fraction 


Scatter  fraction  (SF)  is  defined  as  the  ratio  of  scatter  radiation  to  total  radiation.  Residual 
scatter  fraction  (RSF)  indicates  how  much  of  the  scatter  radiation  remains  after  applying  the 
scatter  compensation  algorithm. 


For  our  imaging  technique,  two  sets  of  images  of  the  phantom  were  obtained:  one  with,  and 
one  without,  a  beam  stop  array.  The  signals  behind  beam  stops  (lead  discs)  comprise  the 
scatter  radiation,  while  the  total  radiation,  which  is  the  sum  of  primary  radiation  and  scatter 
radiation,  will  reach  the  region  without  beam  stops.  We  calculated  the  measured  primary 
radiation  (Pmeasured)  by  subtracting  the  mean  radiation  of  a  region-of-interest  (ROI)  behind  a 
beam  stop  from  the  mean  of  the  same  ROI  location  without  a  beam  stop.  In  the  image 
processed  for  scatter  compensation,  the  mean  of  total  radiation  (T)  in  the  same  ROI  location 
(Testimated)  is  the  sum  of  the  residual  scatter  radiation  and  the  primary  radiation.  Thus, 


RSF 


T  -  P 

estimated  measured 


(7) 


2.3.2  Contrast,  Noise  and  CNR 


Contrast  was  defined  as  the  ratio  of  the  difference  between  the  mean  radiation  value  of  the 
lesion  (Tlesion)  and  the  background  (Tbackground)  to  the  mean  of  the  background,  that  is, 

Contrast  =  T,esi°"  ~Tb“k^'d  ■  (8) 

T 

background 

Noise  was  derived  by  dividing  the  standard  deviation  (STDbackground)  by  the  mean  of  the 
background  radiation  (Thackground): 


STD 

NoiSC  background 

T 

background 

CNR  is  the  ratio  of  the  contrast  to  the  noise,  i.e., 


(9) 


CNR 


Contrast 

Noise 


T  -T 

lesion  background 


STD, 


background 


(10) 


2.3.3  Resolution 

Due  to  the  nonlinearity  of  the  resolution  algorithm,  we  could  not  use  metrics  like  MTF, 
which  are  designed  for  a  linear  system.  Instead,  a  test  bar,  comprised  of  alternating  bright  and 
dark  lines  with  sizes  corresponding  to  Nyquist  frequencies  with  square  wave  function,  was 
embedded  in  the  phantom  image. 

The  contrast  improvement  factor  (CIF),  defined  as  the  ratio  of  the  contrast  after  image 
processing  to  the  initial  contrast,  was  obtained  for  the  test  bar  with  various  initial  contrast 
settings.  A  CIF  of  1  or  greater  was  used  as  the  criterion  for  retaining  the  spatial  resolution.  The 
minimal  initial  contrast  that  the  test  bar  can  allow  with  CIF  of  1  or  greater  was  recorded  as  an 
indication  of  the  effect  of  the  image  processing  on  resolution.  For  the  various  initial  contrasts, 
the  corresponding  CIF  was  computed  arbitrarily  at  iteration  16.  We  determined  the  minimal 
initial  contrast  value  that  has  a  CIF  of  1  or  greater. 

3.  RESULTS 

3.1  Scatter  Compensation  Technique  —  Tissue  Equivalent  Slabs 

Figure  3  and  Table  1  give  the  RSF,  CNR  and  resolution  results  for  MLE  and  MAP  algorithms 
based  on  a  Poisson  noise  model  and  Gaussian  noise  model.  Both  MLE  algorithms  reduced  RSF 
values  to  close  to  zero  and  decreased  CNR  values  from  the  original  unprocessed  value  of  47  to 
slightly  below  40.  The  minimal  contrast  that  is  retainable  during  processing  using  our  Poisson- 
model-based  MAP  algorithm  was  1.8%.  The  MAP  algorithms  were  as  equally  effective  as  their 
MLE  counterparts  in  removing  scattered  radiation  from  the  radiograph;  however,  they  increased 
the  CNR  values  to  56  and  63,  for  the  Poisson  noise  model  and  Gaussian  model,  respectively.  The 
minimal  contrast  retainable  using  the  Gaussian  model  based  MAP  algorithm  was  2.0%,  slightly 
higher  than  that  of  Poisson-model-based  MAP  (1.8%). 
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Figure  3:  Plots  of  (a)  residual  scatter  fraction  and  (b)  contrast  to  noise  ratio  as  the  function  of  iteration 
numbers  between  MLE  and  MAP  estimates  of  scatter  free  image  (bij  based  on  Poisson  and  Gaussian  noise 
models.  While  both  techniques  were  effective  at  removing  scattered  radiation,  the  MAP  based  on  the 
Gaussian  noise  model  showed  greater  CNR  improvement. 


Table  1:  Resolution  results  of  MAP  estimates  based  on  Poisson  noise  model  and  Gaussian  noise  model 
with  the  magnitude  of  scatter  kernel  of  0.52.  The  resolution  results  from  the  two  models  are  close  to  each 
other. 


Poisson  Noise  Model 

Gaussian  Noise  Model 

Minimum  initial  contrast  that  is 
retainable  during  processing 

1.8% 

2.0% 

How  the  magnitude  of  the  scatter  kernel  impacted  the  MAP  algorithm  based  on  the  Gaussian 
noise  model  was  also  investigated.  Figure  4(a)  shows  the  RSF  as  a  function  of  iteration  for 
different  magnitudes.  When  the  magnitude  is  zero,  there  is  no  scatter  reduction  effect.  As  the 
magnitude  increases,  the  steady-state  RSF  decreases.  When  the  magnitude  is  larger  than  the 
measured  value  of  0.52,  the  scattered  radiation  is  overcompensated  such  that  RSF  is  less  than 
zero.  Figure  4(b)  depicts  CNR  results  for  the  same  scatter  kernel  magnitudes.  Overall,  the  larger 
the  magnitude,  the  less  CNR  will  increase.  In  the  case  of  resolution,  the  smaller  the  magnitude, 
the  lower  the  resolution  (Table  2).  For  the  magnitude  of  0.2,  the  minimal  contrast  retainable 
during  processing  was  2.7%.  For  the  magnitude  of  0,  where  there  was  no  scatter  compensation, 
no  contrast  was  retainable  during  processing. 
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Figure  4:  Plots  of  (a)  residual  scatter  fraction  and  (b)  contrast  to  noise  ratio  as  the  function  of  iteration 
numbers  were  shown  for  various  magnitudes  of  scatter  kernel  using  MAP  estimates  of  scatter  free  image 
{ bi j  based  on  the  Gaussian  noise  model.  A  magnitude  of  0.0  corresponds  to  no  scatter  removal,  whereas 
a  magnitude  of  0.65  overcompensates  the  scattered  radiation,  resulting  in  negative  residual  scatter 
fraction  values.  At  each  magnitude  level,  the  contrast  increases  asymptotically.. 


Table  2:  Resolution  results  of  the  Gaussian  noise  model  based  MAP  estimates  with  various  magnitudes  of 
scatter  kernel.  The  larger  the  magnitude  of  scatter  kernel,  the  sharper  the  processed  image  is.  For  the 
magnitude  of  zero,  i.e.,  no  scatter  removal,  the  resolution  is  always  degraded. 


Magnitude 

0.0 

0.2 

0.4 

0.52 

0.65 

Minimal  initial 
contrast 

— 

2.7% 

2.2% 

2.0% 

2.0% 

3.2  Scatter  Compensation  Technique  —  Anthropomorphic  Phantom 

The  scatter  removal  procedure  reduced  SF  of  the  radiograph  acquired  without  a  grid  from  45  %  to 
10%,  the  level  that  an  anti-scatter  grid  achieves  (Figure  5,  Table  3).  At  the  same  time,  the 
procedure  improved  the  CNR  to  around  twice  the  value  on  the  image  acquired  with  a  grid. 


(a)  with  grid  (b)  without  grid  (c)  without  grid;  scatter  reduction 


Figure  5:  Radiographs  of  the  breast  anthropomorphic  phantom,  (a)  with  an  anti-scatter  grid,  (b)  without 
an  anti-scatter  grid,  and  (c)  without  an  anti-scatter  grid  and  with  scatter  reduction. 


Table  3:  Corresponding  residual  scatter  fraction  and  contrast  to  noise  ratio  results  for  the  three  images 
shown  in  Figure  5. 


With  grid 

w/o  grid 

w/o  grid;  scatter  reduction 

RSF 

11% 

45% 

10% 

CNR 

7.04 

6.99 

15.29 

4.  DISCUSSION 

In  this  study,  both  MLE  and  MAP  estimates  of  the  scatter  free  image  were  derived  based  on  a 
novel  Gaussian  noise  model  for  energy-integrating  detectors.  The  preliminary  results  were 
obtained  on  the  two  types  of  phantoms  (tissue  equivalent  slabs  and  a  breast  anthropomorphic 
phantom)  obtained  on  a  full-field  digital  mammography  system.  Both  MLE  and  MAP  techniques 
were  effective  in  removing  the  scattered  radiation,  though  MAP  outperformed  MLE  in  CNR.  For 
the  specific  phantom  and  imaging  condition,  the  MAP  of  the  Gaussian  noise  model  outperformed 
the  MAP  of  the  Poisson  noise  model. 

The  next  phase  of  this  research  will  include  a  comprehensive  evaluation  of  the  scatter  reduction 
technique  on  images  acquired  on  a  FFDM  system.  Also  the  technique  can  be  applied  on  dedicated 
breast  CT  data  for  scatter  reduction. 
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6.  APPENDIX 

Derivation  of  MLE  Algorithm  based  on  Gaussian  Noise  Model 


Due  to  the  convolution  operation,  the  estimation  of  B  =  {bj;  i=l...N}  directly  from  Y  does  not 
have  a  simple  analytic  form.  The  MLE  of  B  is  thus  derived  through  the  EM  algorithm  as  follows. 

Treat  the  measured  Y  =  {y;,  i=l,...,N}  as  an  incomplete  dataset,  and  unobserved  (D,S)  =  {(d;,Si), 
i=l,...,N}  as  a  complete  dataset.  The  d;  ’s  and  s;  ’s  given  B  are  mutually  independent,  therefore 
the  complete  data  likelihood  is 


pc(D,S\BAoa2,oi22-i  =  l,...,N))  =  Y[ 


-(dj-bj)  I Icfji 


-{Sj-(b**p)j)  12a j2 


h  J2jta, 


]2 KO 


(Al) 


Assuming  {on2,  cfo2 ;  i=l,...,N}  are  known,  we  can  obtain  the  complete  data  log  likelihood  by 
taking  the  logarithm  on  both  sides. 


(B\D,S)  =  J/[~(dJ  -bi)2/2o]l2-(Sj-(b**p)j)2/2oii  -  log  ^2aa/  -  log  ^2jio  J22  ]  • 


(A2) 


The  EM  algorithm  is  comprised  of  two  steps:  one,  the  E-step,  where  the  expectation  of  the 
complete  data  log  likelihood  with  respect  to  the  present  estimate  of  B  is  computed,  and  two,  the 
M-step,  where  a  new  estimate  of  B  is  obtained  which  will  maximize  the  computed  expectation  in 
the  E-step. 

First,  employ  the  E-step: 


Q(B\Bm) 

=  E[Lc(B\D,S)\Y,B{n)]  (A3) 

=  {-(£>/  -2d^bj)l2o  2  ~[(B**P) j)2  -2sf")(B**P)j]/2oj22 

+  terms  •  independent  •  of  •  B), 

wher  edjM=E[dj\Y,B^]  (A4) 

s/n)  =E[Sj  I  Y,BW] 

Second,  the  M-step  to  find  B(n+I)  that  will  maximize  Q(BIB(n)): 


dQ(B  \Bn  )  9  ,  (n)s/9  2  Vn^R**pT 

=  0  =  ~(2bk  -  2dk  )/2 okl  ~2J[2(B**P)jpjk- 


db,. 


2spPjk]/2oj22- 


7=1 


Solving  the  above  equation  for  bk  gives 


bp"  =  dkM -akl2yPjk  ■[(B(n+l)**P)j  -spy 


Uj2  . 


(A5) 


Using  B(n)  to  approximate  BU1+1J  in  the  right-hand  side  yields 


,(n+l) 


bp"-dkM  _au2yPjk  .[(/}'">  **P)j -SjM]/a„2. 


(A6) 


As  a  good  estimate  of  the  primary  image  is  formed,  (B<n>  **P) .  -  s  M  =  0 ,  then 


(A7) 

The  same  apparent  form  was  obtained  for  Poisson  noise  model  in  Reference  2.  But  due  to  the 
different  statistical  models,  the  actual  forms  of  dk(n)  are  different  and  so  is  the  iterative  formula  for 
bk. 

Equation  (A7)  combines  with  equation  (A4)  to  give  the  following  updated  equation: 

bp"  =  bp'  +  wk  •  [yk-(bp  +  (B,">**P)k)\,  (A8) 


where 


(A9) 
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